Flowcytometer Biomass Calculation

Introduction

Details the calculation of the biomass and the problem occuring with FSC.A == 0

Code

library(LEEF.analysis)
library(LEEF.measurement.flowcytometer)
library(scattermore)

Calculate measurements

Code

if (!dir.exists(file.path(params$output_dir, "pre_processed"))){
  dir.create(file.path(params$output_dir, "pre_processed"), showWarnings = FALSE, recursive = TRUE)
  
  file.copy(
    file.path(params$pre_processed_dir, paste0("LEEF.fast.flowcytometer.", as.character(params$timestamp))),
    file.path(params$output_dir, "pre_processed"),
    recursive = TRUE,
    copy.mode = FALSE
  )
  
  file.rename(
    from = file.path(params$output, "pre_processed", paste0("LEEF.fast.flowcytometer.", as.character(params$timestamp))),
    to   = file.path(params$output, "pre_processed", "flowcytometer")
  )
}

if (!dir.exists(file.path(params$output_dir, "extracted_raw"))){
  dir.create(file.path(params$output_dir, "extracted_raw"), showWarnings = FALSE, recursive = TRUE)
  
  LEEF.measurement.flowcytometer::extractor_flowcytometer(
    input = file.path(params$output, "pre_processed"),
    output = file.path(params$output_dir, "extracted_raw"),
    raw = TRUE
  )
}

if (!dir.exists(file.path(params$output_dir, "extracted_log10"))){
  dir.create(file.path(params$output_dir, "extracted_log10"), showWarnings = FALSE, recursive = TRUE)
  
  LEEF.measurement.flowcytometer::extractor_flowcytometer(
    input = file.path(params$output, "pre_processed"),
    output = file.path(params$output_dir, "extracted_log10"),
    raw = FALSE
  )
}

Read Data

Size Standard

Code

ss_replacement <- list()
ss_replacement$data <- data.frame(
  diameter_micrometer = c(1L, 3L, 10L), 
  mean_FSC.A = c(96374L, 653207L, 1855322L), 
  mean_FSC.H = c(100606L, 661955L, 1252168)
)

ss_replacement$fit.A <- lm(ss_replacement$data$diameter_micrometer ~ ss_replacement$data$mean_FSC.A)
ss_replacement$slope_A <- ss_replacement$fit.A$coefficients[[2]]
ss_replacement$intercept_A <- ss_replacement$fit.A$coefficients[[1]]

ss_replacement$fit.H <- lm(ss_replacement$data$diameter_micrometer ~ ss_replacement$data$mean_FSC.H)
ss_replacement$slope_H <- ss_replacement$fit.H$coefficients[[2]]
ss_replacement$intercept_H <- ss_replacement$fit.H$coefficients[[1]]


ss_repaired <- list()
ss_repaired$data <- data.frame(
  diameter_micrometer = c(1L, 3L, 10L), 
  mean_FSC.A = c(102045L, 641752L, 1959230L), 
  mean_FSC.H = c(106576L, 648403L, 1303073L)
)

ss_repaired$fit.A <- lm(ss_repaired$data$diameter_micrometer ~ ss_repaired$data$mean_FSC.A)
ss_repaired$slope_A <- ss_repaired$fit.A$coefficients[[2]]
ss_repaired$intercept_A <- ss_repaired$fit.A$coefficients[[1]]

ss_repaired$fit.H <- lm(ss_repaired$data$diameter_micrometer ~ ss_repaired$data$mean_FSC.H)
ss_repaired$slope_H <- ss_repaired$fit.H$coefficients[[2]]
ss_repaired$intercept_H <- ss_repaired$fit.H$coefficients[[1]]

Log10 transformed

Code

datadir <- file.path(params$output_dir, "extracted_log10", "flowcytometer")

write.csv(
  data.frame(
    bacteria1 = c(3.477121255, 3.477121255, 7.204119983, 7.204119983),
    bacteria2 = c(0,           3.301029996, 6.447158031, 0),
    LNA = c(3,4.579784, NA, NA),
    MNA = c(4.579785, 4.954243, NA, NA),
    HNA = c(4.954244, 7, NA, NA),
    algae1 = c(3, 3, 7, 7),
    algae2 = c(4, 7, 7, 4)
  ),
  row.names = FALSE,
  file = file.path(datadir, "gates_coordinates.csv")
)

traits_log10 <- LEEF.measurement.flowcytometer::extract_traits(
  input = datadir,
  particles = params$particles,
  metadata_flowcytometer = read.csv(file.path(datadir, "metadata_flowcytometer.csv"))
)

{"timestamp": "2023-03-07T13:09:54+0100", "log_lvl": "INFO", "log_msg": "   reading data  ..."}

   reading data  ...

{"timestamp": "2023-03-07T13:09:56+0100", "log_lvl": "INFO", "log_msg": "   gating bacteria ..."}

   gating bacteria ...

{"timestamp": "2023-03-07T13:10:11+0100", "log_lvl": "INFO", "log_msg": "   done"}

   done

{"timestamp": "2023-03-07T13:10:11+0100", "log_lvl": "INFO", "log_msg": "########################################################"}

########################################################

Code

paste0("Count       : ",  traits_log10[[1]]$FSC.A |> length())

[1] "Count       : 3780167"

Code

paste0("Range       : ",  paste0(traits_log10[[1]]$FSC.A |> range(), collapse = " -- "))

[1] "Range       : 0 -- 7.22471987004958"

Code

paste0("Portion == 0: ", ((traits_log10[[1]]$FSC.A == 0) |> sum()) / (traits_log10[[1]]$FSC.A |> length()))

[1] "Portion == 0: 0.120381718585449"

Code

traits_log10[[1]]$FSC.A |> table() |> head()

                0 0.301029995663981 0.477121254719662 0.602059991327962 
           455063               946               942              1863 
0.698970004336019 0.778151250383644 
              954               946 

Raw

Code

datadir <- file.path(params$output_dir, "extracted_raw", "flowcytometer")

write.csv(
  data.frame(
    bacteria1 = 10^c(3.477121255, 3.477121255, 7.204119983, 7.204119983),
    bacteria2 = 10^c(0,           3.301029996, 6.447158031, 0),
    LNA = c(3,4.579784, NA, NA),
    MNA = c(4.579785, 4.954243, NA, NA),
    HNA = c(4.954244, 7, NA, NA),
    algae1 = c(3, 3, 7, 7),
    algae2 = c(4, 7, 7, 4)
  ),
  row.names = FALSE,
  file = file.path(datadir, "gates_coordinates.csv")
)


traits_raw <- LEEF.measurement.flowcytometer::extract_traits(
  input = datadir,
  particles = params$particles,
  metadata_flowcytometer = read.csv(file.path(datadir, "metadata_flowcytometer.csv"))
)

{"timestamp": "2023-03-07T13:10:13+0100", "log_lvl": "INFO", "log_msg": "   reading data  ..."}

   reading data  ...

{"timestamp": "2023-03-07T13:10:15+0100", "log_lvl": "INFO", "log_msg": "   gating bacteria ..."}

   gating bacteria ...

{"timestamp": "2023-03-07T13:10:26+0100", "log_lvl": "INFO", "log_msg": "   done"}

   done

{"timestamp": "2023-03-07T13:10:26+0100", "log_lvl": "INFO", "log_msg": "########################################################"}

########################################################

Code

paste0("Count       : ",  traits_raw[[1]]$FSC.A |> length())

[1] "Count       : 3712417"

Code

paste0("Range       : ",  paste0(traits_raw[[1]]$FSC.A |> range(), collapse = " -- "))

[1] "Range       : 0 -- 16777215"

Code

paste0("Portion == 0: ", ((traits_raw[[1]]$FSC.A == 0) |> sum()) / (traits_raw[[1]]$FSC.A |> length()))

[1] "Portion == 0: 0.120776033511322"

Code

traits_raw[[1]]$FSC.A |> table() |> head()

     0      1      2      3      4      5 
448371    933    936    934   1838    944 

Some Plots

Code

i0 <- traits_raw[[1]]$FSC.A == 0
i1 <- traits_raw[[1]]$FSC.A == 1

Size Standards

Repaired

Code

plot(
  ss_repaired$data$mean_FSC.A, 
  ss_repaired$data$diameter_micrometer, 
  xlim = c(0, 3000000), 
  ylim = c(0, 20),
  main = "Repaired - FSC.A"
)
abline(a = ss_repaired$intercept_A, b = ss_repaired$slope_A)
abline(v = min(traits_raw[[1]]$FSC.A))
abline(v = max(traits_raw[[1]]$FSC.A))
abline(h = 0)

Code

plot(
  ss_repaired$data$mean_FSC.H, 
  ss_repaired$data$diameter_micrometer, 
  xlim = c(0, 1500000), 
  ylim = c(0, 20),
  main = "Repaired - FSC.H"
)
abline(a = ss_repaired$intercept_H, b = ss_repaired$slope_H)
abline(v = min(traits_raw[[1]]$FSC.H))
abline(v = max(traits_raw[[1]]$FSC.H))
abline(h = 0)

Replacement

Code

plot(
  ss_replacement$data$mean_FSC.A, 
  ss_replacement$data$diameter_micrometer, 
  xlim = c(0, 3000000),
  ylim = c(0, 20),
  main = "Replacement - FSC.A"
)
abline(a = ss_replacement$intercept_A, b = ss_replacement$slope_A)
abline(v = min(traits_raw[[1]]$FSC.A))
abline(v = max(traits_raw[[1]]$FSC.A))
abline(h = 0)

Code

plot(
  ss_replacement$data$mean_FSC.H, 
  ss_replacement$data$diameter_micrometer, 
  xlim = c(0, 1500000), 
  ylim = c(0, 20),
  main = "Replacement - FSC.H"
)
abline(a = ss_replacement$intercept_H, b = ss_replacement$slope_H)
abline(v = min(traits_raw[[1]]$FSC.H))
abline(v = max(traits_raw[[1]]$FSC.H))
abline(h = 0)

Scatterplots

FSC-A auf der x-Achse und SSC-A auf der y-Achse

Code

scattermoreplot(traits_raw[[1]]$FSC.A, traits_raw[[1]]$SSC.A, xlab = "FSC.A", ylab = "SSC.A")

Code

x <- traits_log10[[1]]$FSC.H
x[x==0] <- 1
x <- log10(x)
scattermoreplot(traits_log10[[1]]$FSC.A, x, xlab = "log10(FSC.A)", ylab = "log10(SSC.A)")

FSC-A auf der x-Achse und FL1-A auf der y-Achse

Code

scattermoreplot(traits_raw[[1]]$FSC.A, traits_raw[[1]]$FL1.A, xlab = "FSC.A", ylab = "FL1.A")

Code

x <- traits_log10[[1]]$FL1.A
x[x==0] <- 1
x <- log10(x)
scattermoreplot(traits_log10[[1]]$FSC.A, x, xlab = "log10(FSC.A)", ylab = "log10(FL1.A)")

FSC-A auf der x-Achse und FL3-A auf der y-Achse

Code

scattermoreplot(traits_raw[[1]]$FSC.A, traits_raw[[1]]$FL3.A, xlab = "FSC.A", ylab = "FL3.A")

Code

x <- traits_log10[[1]]$FL3.A
x[x==0] <- 1
x <- log10(x)
scattermoreplot(traits_log10[[1]]$FSC.A, x, xlab = "log10(FSC.A)", ylab = "log10(FL3.A)")

FSC-A auf der x-Achse und FSC-H auf der y-Achse

Code

scattermoreplot(traits_raw[[1]]$FSC.A, traits_raw[[1]]$FSC.H, xlab = "FSC.A", ylab = "FSC.H")

Code

x <- traits_log10[[1]]$FSC.H
x[x==0] <- 1
x <- log10(x)
scattermoreplot(traits_log10[[1]]$FSC.A, x, xlab = "log10(FSC.A)", ylab = "log10(FSC.H)")

And some more density plots

Raw FSC.A and FSC.H

Code

traits_raw[[1]]$FSC.A |> density(bw = 10, to = 5000) |> plot(main = "FSC.A")

Code

traits_raw[[1]]$FSC.H |> density(bw = 10, to = 10000) |> plot(main = "FSC.H")

log10 FSC.A and FSC.H

Code

traits_log10[[1]]$FSC.A |> density(bw = 0.02) |> plot(main = "log10(FSC.A)")

Code

x <- traits_log10[[1]]$FSC.H
x[x==0] <- 1
x |> log10() |> density(bw = 0.02) |> plot(main = "log10(FSC.H)")

• BLACK Line: All particles with FSC.A == 0 (will tbe transformed to 0 using log10)
• RED Line: All particles with FSC.A == 1 (will tbe transformed to 0 using log10)
• Black Line: All particles with FSC.A != 0

SSC.A

Code

traits_raw[[1]]$SSC.A[i0] |> density(to = 10000, bw = 1) |> plot()
traits_raw[[1]]$SSC.A[i1] |> density(to = 10000, bw = 1) |> lines(col = "red")
traits_raw[[1]]$SSC.A[!i0] |> density(to = 10000, bw = 1) |> lines(col = "green")

FL1.A

Code

traits_raw[[1]]$FL1.A[i0] |> density(to = 10000, bw = 1) |> plot()
traits_raw[[1]]$FL1.A[i1] |> density(to = 10000, bw = 1) |> lines(col = "red")
traits_raw[[1]]$FL1.A[!i0] |> density(to = 10000, bw = 1) |> lines(col = "green")

FL2.A

Code

traits_raw[[1]]$FL2.A[i0] |> density(to = 10000, bw = 1) |> plot()
traits_raw[[1]]$FL2.A[i1] |> density(to = 10000, bw = 1) |> lines(col = "red")
traits_raw[[1]]$FL2.A[!i0] |> density(to = 10000, bw = 1) |> lines(col = "green")

FL3.A

Code

traits_raw[[1]]$FL3.A[i0] |> density(to = 10000, bw = 1) |> plot()
traits_raw[[1]]$FL3.A[i1] |> density(to = 10000, bw = 1) |> lines(col = "red")
traits_raw[[1]]$FL3.A[!i0] |> density(to = 10000, bw = 1) |> lines(col = "green")

FL4.A

Code

traits_raw[[1]]$FL4.A[i0] |> density(to = 2000, bw = 1) |> plot()
traits_raw[[1]]$FL4.A[i1] |> density(to = 2000, bw = 1) |> lines(col = "red")
traits_raw[[1]]$FL4.A[!i0] |> density(to = 2000, bw = 1) |> lines(col = "green")

FL2.A

Code

traits_raw[[1]]$FL2.A[i0] |> density(to = 10000, bw = 1) |> plot()
traits_raw[[1]]$FL2.A[i1] |> density(to = 10000, bw = 1) |> lines(col = "red")
traits_raw[[1]]$FL2.A[!i0] |> density(to = 10000, bw = 1) |> lines(col = "green")

FSC.H

Code

traits_raw[[1]]$FSC.H[i0] |> density(to = 10000, bw = 1) |> plot()
traits_raw[[1]]$FSC.H[i1] |> density(to = 10000, bw = 1) |> lines(col = "red")
traits_raw[[1]]$FSC.H[!i0] |> density(to = 10000, bw = 1) |> lines(col = "green")

SSC.H

Code

traits_raw[[1]]$SSC.H[i0] |> density(to = 10000, bw = 1) |> plot()
traits_raw[[1]]$SSC.H[i1] |> density(to = 10000, bw = 1) |> lines(col = "red")
traits_raw[[1]]$SSC.H[!i0] |> density(to = 10000, bw = 1) |> lines(col = "green")

FL1.H

Code

traits_raw[[1]]$FL1.H[i0] |> density(to = 10000, bw = 1) |> plot()
traits_raw[[1]]$FL1.H[i1] |> density(to = 10000, bw = 1) |> lines(col = "red")
traits_raw[[1]]$FL1.H[!i0] |> density(to = 10000, bw = 1) |> lines(col = "green")

FL2.H

Code

traits_raw[[1]]$FL2.H[i0] |> density(to = 10000, bw = 1) |> plot()
traits_raw[[1]]$FL2.H[i1] |> density(to = 10000, bw = 1) |> lines(col = "red")
traits_raw[[1]]$FL2.H[!i0] |> density(to = 10000, bw = 1) |> lines(col = "green")

FL3.H

Code

traits_raw[[1]]$FL3.H[i0] |> density(to = 10000, bw = 1) |> plot()
traits_raw[[1]]$FL3.H[i1] |> density(to = 10000, bw = 1) |> lines(col = "red")
traits_raw[[1]]$FL3.H[!i0] |> density(to = 10000, bw = 1) |> lines(col = "green")

FL4.H

Code

traits_raw[[1]]$FL4.H[i0] |> density(to = 10000, bw = 1) |> plot()
traits_raw[[1]]$FL4.H[i1] |> density(to = 10000, bw = 1) |> lines(col = "red")
traits_raw[[1]]$FL4.H[!i0] |> density(to = 10000, bw = 1) |> lines(col = "green")

Width

Code

traits_raw[[1]]$Width[i0] |> density(to = 500, bw = 1) |> plot()
traits_raw[[1]]$Width[i1] |> density(to = 500, bw = 1) |> lines(col = "red")
traits_raw[[1]]$Width[!i0] |> density(to = 500, bw = 1) |> lines(col = "green")

Now to the sizes

FSC.A based

Code

#|

i <- traits_log10[[1]]$FSC.A > 1

traits_log10[[1]]$length_A[i] <- 10^(traits_log10[[1]]$FSC.A[i]) * ss_repaired$slope_A + ss_repaired$intercept_A
traits_log10[[1]]$volume_A <- 4/9 * pi * traits_log10[[1]]$length_A^3

traits_log10[[1]]$length_A |> density(na.rm = TRUE) |> plot(main = "FSC.A - length in micro meter")

Code

traits_log10[[1]]$length_A |> density(na.rm = TRUE, to = 0.5) |> plot(main = "FSC.A - length in micro meter (l < 0.5")

Code

traits_log10[[1]]$length_A |> density(na.rm = TRUE, to = 0.3) |> plot(main = "FSC.A - length in micro meter (l < 0.2)")

FSC.H based

Code

#|

i <- traits_log10[[1]]$FSC.H != 0

traits_log10[[1]]$length_H <- traits_log10[[1]]$FSC.H * ss_repaired$slope_H + ss_repaired$intercept_H
traits_log10[[1]]$volume_H <- 4/9 * pi * traits_log10[[1]]$length_H^3

traits_log10[[1]]$length_H |> density(na.rm = TRUE) |> plot(main = "FSC.H - length in micro meter")

Code

traits_log10[[1]]$length_H |> density(na.rm = TRUE, to = 2) |> plot(main = "FSC.H - length in micro meter (l < 2)")

Code

traits_log10[[1]]$length_H |> density(na.rm = TRUE, to = 0.5) |> plot(main = "FSC.H - length in micro meter (l < 0.5)")

Code

(traits_log10[[1]]$length_H - traits_log10[[1]]$length_A) |>
  density(na.rm = TRUE) |>
  plot(main = "length_H - length_A")

Code

(traits_log10[[1]]$length_H - traits_log10[[1]]$length_A) |>
  density(na.rm = TRUE, from = -1, to = 1) |>
  plot(main = "length_H - length_A")